Competing many-body instabilities and unconventional superconductivity in graphene 
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The band structure of graphene exhibits van Hove singularities (VHS) at doping x = ±1/8 
away from the Dirac point. Near the VHS, interactions effects, enhanced due to the large density 
of states, can give rise to various many-body phases at experimentally accessible temperatures. 
We study the competition between different many-body instabilities in graphene using functional 
renormalization group (FRG). We predict a rich phase diagram, which, depending on long range 
hopping as well as screening strength and absolute scale of the Coulomb interaction, contains a 
d + id-wave superconducting (SC) phase, or a spin density wave phase at the VHS. The d + id state 
is expected to exhibit quantized charge and spin Hall response, as well as Majorana modes bound to 
vortices. In the vicinity of the VHS, we find singlet d-|- id-wave as well as triplet /-wave SC phases. 

PACS numbers: 73.22.Pr,74.70.Wz,74.20.Mn 



Introduction. Graphene, a monolayer of carbon, hosts 
a two-dimensional electronic system (2DES) with unique 
properties P]. In particular, the Coulomb interaction 
plays an important role in graphene [2], giving rise to 
interesting many-body phenomena, including marginal 
Fermi- liquid behavior |3] , energy-dependent renormaliza- 
tion of the Fermi velocity [3], as well many-body states 
in the quantum Hall effect regime [I] . 

Experimentally, graphene offers a high degree of tun- 
ability. Most importantly, carrier density can be con- 
trolled in a broad range. Near the Dirac point (doping 
level of electrons x = 1/2), such control is achieved by 
backgates and local top gates [1] . Recently, it was demon- 
strated that chemical doping [5] and electrolytic gat- 
ing ^ enable doping graphene far away from the Dirac 
point, where the band structure is no longer Dirac-like. 
In particular, the density can be tuned to the vicinity of 
the van Hove singularities (VHS) in the band structure, 
which occur at doping values x = 3/8, 5/8. In the case of 
chemical doping, the dopants form a superlattice on top 
of the graphene sheet. This strongly reduces the amount 
of disorder induced by doping. Furthermore, the spac- 
ing of the superlattice is so large that hybridization in 
the dopant layer can be neglected, and hence transport 
measurements of the graphene sample remain unaffected. 

Before the strong doping of graphene has been re- 
cently accomplished experimentally, superconductivity 
had been predominantly studied around the Dirac point, 
including p + ip-wsve from electron-phonon or plas- 
monic [7 as well as / or d -I- id-wave from electron- 
electron |H] interaction effects. Mean-field treatments [51 
110] were found to be unreliable as they arrive at unrealis- 
tic Tc > 1000 K, with only slightly better results for vari- 
ational approaches [TT]. SC has not been observed exper- 
imentally in this regime, due to small electronic density 
of states as well as weak phonon effects |12) . 

Near a VHS, opposed to the Dirac point regime, elec- 
tronic interaction effects are expected to be strongly en- 



hanced due to the logarithmically diverging density of 
states and near-nested Fermi-surface [5|. In this regime, 
many-body states with appreciable critical temperatures 
may arise. Possible candidate states include charge- 
density wave (CDW), a spin-density wave (SDW), or 
a SC state. Generally, a subtle interplay of kinetic 
and interaction parameters is expected to decide which 
many body instability is preferred at the VHS. For 
graphene, the additional complication arises that as the 
band width (^ 17eV) is of the order of the interaction 
scale (~ lOeV), graphene cannot be suitably described 
from the viewpoint of strict weak coupling approaches, 
and adopting a picture of intermediate coupling is nec- 
essary. Rephrased in terms of diagrammatic expansions 
starting from the non-interacting problem, this amounts 
to investigating the importance of leading and sublead- 
ing divergent classes of diagrams. In particular, this is 
relevant for the competition between magnetic and SC 
phases in this kind of systems, one recent example of 
which have been the iron pnictides [HlITi] . 

Main results. In this Letter, we use the functional 
renormalization group (FRG) method [T5HT5] to study 
the competition between many-body states in graphene 
doped to the vicinity of the VHS, and attempt to an- 
alyze this problem at a level which provides a detailed 
connection to the experimental setup. From our analy- 
sis we obtain a rich phase diagram which, depending on 
the range of chosen kinetic and interaction parameters, 
contains magnetic and different SC phases, summarized 
in Fig. [l] For a certain range of parameters, we find 
a. d + id SC phase which has been previously studied 
by RPA [51 [12], and, very recently, perturbative three- 
patch renormalization group (3RG) [2D]. To analyze all 
possible many-body phases and their dependence on the 
system parameters, FRG provides a systematic unbiased 
summation of diagrams in both particle-particle channels 
and particle-hole channels as well as vertex corrections, 
and keeps track of the whole Fermi surface [Fig. [2^]. We 
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FIG. 1. (Color online). Schematic phase diagram displaying 
the critical instability scale ~ Tc as a function of dop- 
ing. At the van Hove singularity (VHS, light shaded (orange) 
area), d + id competes with spin density wave (SDW) (left 
flow picture: dominant d + id instability for Uo = lOeV and 
the band structure in [5]). Away from the VHS (dark shaded 
(blue) area), Ac drops and whether the d+id or /-wave SC in- 
stability is preferred depends on the long-rangedness of inter- 
action (right flow picture: Ui/Uq — 0.45 and U2/U0 ~ 0.15). 




FIG. 2. (Color online), (a) Band structure of graphene once 
for h = 2.8eV (red) and ti = 2.8, t2 = 0.7, = 0.02eV 
(black), (b) Brillouin zone displaying the Fermi surface near 
the van Hove point (dashed blue level in (a), 96 patches used 
in the FRG and the nesting vector, and the partial nesting 
vectors, (c) Density of states for both band structures in (a). 
The inset show the position shift of Fermi surface nesting 
(dashed vertical lines) versus the VHS peak. 



investigate in detail how different band structure param- 
eters affect the phase diagram. We find that rather smaU 
variations of the longer range hopping parameters such as 
next nearest (12) and next-next-nearest (^3) hopping can 
shift the position of perfect Fermi surface nesting against 
the VHS [Fig. [2], which significantly influences the com- 
petition between magnetism and SC. Moreover, in par- 
ticular away from the exact VHS, the reduced screening 
of the Coulomb interaction does not justify the assump- 
tion of a local Hubbard model description. For this case, 
we find that only a small fraction of longer-ranged Hub- 
bard interaction [3T] can significantly change the phase 
diagram, as CDW fluctuations become more competitive 
to SDW fluctuations, and a triplet SC phase can appear. 
In particular, we study how the Cooper pairing in the 
different SC phases responds to differently long-ranged 
Hubbard interactions. Our results suggest that in ex- 
periment, modiflcations of the band structure as well as 
changing the dielectric environment of the graphene sam- 
ple would enable the realization of different many-body 
states and possible phase transitions between them. 

Model. We consider the vr band structure of graphene 
approximated by a tight binding model including up to 
3rd nearest neighbors on the hexagonal lattice: 

Ho ^[tiJ2Yl ^U^J:'^ + ^2 XI H ^US.'T 

{^■J) ((«J'>> 

where n = Y,^^^ Ui^^ = ^ cl^c^ ^, and cj denotes the 
electron annihilation operator of spin a =t,4- at site i. 



The resulting band structure is a two band model due 
to two atoms per unit cell [Fig. [2]. There are certain 
uncertainties about the most appropriate tight binding flt 
for graphene, in particular as it concerns the longer range 
hybridization integrals [1] [21] . For dominant ti , the band 
structure features a van Hove singularity (VHS) at a; = 
3/8,5/8. Constraining ourselves to the electron-doped 
case, the x = 5/8 electron-like Fermi surface is shown 
in Fig. [2]d. As depicted, this is the regime of largely 
enhanced density of states which we investigate in the 
following. For t2 = ^3 = [red curve in Fig. [2], the VHS 
coincides with the partial nesting of different sections of 
the Fermi surface for Q = (0, 27r/\/3), (tt, 7r/-\/3), and 
(tt, — 7r/%/3). For a realistic band structure estimate with 
flnite t2 and 5 [black curve in Fig. [2], this gives a 
relevant shift of the perfect nesting position versus the 
VHS as well as density of states at the VHS, and affects 
the many-body phase found there. 

We assume Coulomb interactions represented by a long 
range Hubbard Hamiltonian '21 

+ \u2 X ?^^,<TnJ>', (2) 

<T,(j' 

where ?7o...2 parametrizes the Coulomb repulsion scale 
from onsite to the second nearest neighbor interaction. 
It depends on the density of states how strongly the 
Coulomb interaction is screened. At the VHS, we as- 
sume perfect screening and consider Uq only, while away 
from the VHS, we investigate the phenomenology of tak- 
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ing Ui and U2 into consideration. The typical scale of 
the effective Uq has been found to be ~ lOeV < W [21], 
where W ~ 17eV is the kinetic bandwidth. 

Method. We employ the FRG and study how the renor- 
malized interaction described by the 4-point function 
(4PF) evolves under integrating high energy fermionic 
modes: VX(A;i; fc2; fea; k4)cl^^cl^gC^^^c^^~, where the flow 
parameter is the IR cutoff A approaching the Fermi sur- 
face, and with fei to the incoming and outgoing mo- 
menta. Within the numerical treatment, the fc's are dis- 
cretized to take on the values representing the different 
patches of the Brillouin zone. Fig.jSja shows a 96 patching 
scheme. We checked for selected representative scenar- 
ios that our results are converged against supercomputer 
simulations with 192 patch resolution. The starting con- 
ditions of the RG are given by the bandwidth W serving 
as an UV cutoff, with the bare initial interactions for 
the 4PF. Due to the spin rotational invariance of inter- 
actions (we neglect spin-orbit coupling in our analysis), 
we constrain ourselves to the 5*^ = subspace of incom- 
ing momenta fci,fc2 (and outgoing k^jki) and generate 
the singlet and triplet channel by symmetrization and 
antisymmetrization of the 4PF Va [T7j. The diverging 
channels of the 4PF under the flow to the Fermi surface 
signal the nature of the instability, and the corresponding 
Ac, as a function of some given system parameter such 
as doping, gives the same qualitative behavior as Tc- At 
a cutoff scale where the leading instability starts to di- 
verge, we then decompose the different channels such as 
SC or SDW into different eigenmode contributions and 
obtain the form factors associated with the different in- 
stabilities [T31I23]. 

Phase diagram. The phase diagram as a function of 
doping, obtained for realistic microscopic kinetic and in- 
teraction parameters [H] is shown in Fig. [l] At the 
VHS [orange-shaded area in Fig. [l] , the density of states 
is so large that a local Hubbard description is appropri- 
ate. For Uq ~ lOeV, we find that the d + id SC insta- 
bility is dominant, assuming finite hopping parameters 
t2 = 0.7eV and ^3 = 0.02eV [5]. (The result is rather 
similar for the values of [32].) Only at scales of Uq > 18 
eV the SDW becomes dominant for this scenario. Note, 
however, that only small modifications of the band struc- 
ture can strongly affect the competition of SDW and SC 
at the VHS: when t2 is reduced, the system gets more bi- 
ased to the SDW, as the SDW fluctuations in the nesting 
channel get enhanced. For ti only [red curve in Fig. [2^], 
the SDW already wins for Uq > 8.5 eV. As we move away 
from the VHS [blue-shaded area in Fig. [l] , details of the 
band structure become less relevant, and we note that 
the critical instability scale Ac drops stronger towards the 
Dirac point than away from it, mainly as a consequence 
of the reduced density of states. As SDW fluctuations 
are weakened, SC phases become dominant. Still assum- 
ing rather local Coulomb interactions {Ui/Uq < 0.40), we 
find that the system still favors the d+id SC state. AUow- 




FIG. 3. (Color online). d^2_y2 and d^y-wave solutions for (a) 
Uo = lOeV only (and (b) Ui/Uo = 0.4. (al), (a2) and (bl), 
(b2) show the form factors of d^2_y2 and d^y plotted along the 
Fermi surface according to patch indices defined in Fig.|2]D, as 
well as the real space pair amplitude patterns. The solutions 
change from (a) to (b). The analytic form factors given in 
the text (red) fit the numerical data (black). (a3) and (b3) 
show the gap profile of d -f id along the Fermi surface (actual 
connection to experimental energy scale can still vary by a 
global factor). The gap anisotropy increases from (a) to (b). 



ing for more long-ranged Hubbard interactions, however, 
the picture changes: the CDW fluctuations are compara- 
ble to the SDW fluctuations which would bias the system 
towards singlet SC, and triplet /-wave pairing becomes 
competitive. 

d + id-wave phase. Let us analyze the d-wave SC 
phase at the VHS {Uq only) in more detail. The honey- 
comb lattice is characterized by Cgu symmetry about the 
center of hexagons, and the SC order parameter trans- 
forms as one of the irreducible representations. dx2-y2 
and fia;j,-wave follow the two-dimensional E2 representa- 
tion and are hence degenerate. The different form fac- 
tors are plotted in Fig. [3^. We find that the numer- 
ical solutions can be fit to f[dx2_y2] = 2cos{^/?>ky) — 
cos[(-\/3fcy — 2>kx)/2] — cos[{y/iky -t- 'ikx)/'2\ and f[dxy\ = 
cos[('\/3/cy — ^kx)/2] — cos[(\/3fcj, -I- ikx)/2\. From the 
Fourier transform of the momentum-resolved form fac- 
tors /(fe) along the Fermi surface we also obtain the pair- 
ing amplitudes of the real space SC pairing function [23] 
[Fig. [3] . The Cooper pairing emerges on nearest neigh- 
bors of the same hexagonal sublattice. As we move to 
the broader vicinity of the VHS where we assume longer 
range Hubbard interaction, the form factors retain the d- 
wave E2 representation, while the Cooper pair wave func- 
tion changes as shown in Fig.[3}D {Ui/Uq = 0.4, U2/U0 = 
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0.25). There, the form factor fits change to f[dx2-y2] = 
2cos{3kx) — cos[(3-\/3/cy — 3fcx)/2] — cosKSVSky + Skx)/^] 
and f[dxy] = cos[{3V3ky--3kx) /2]—cos[{3VSky+3kx) /2], 
corresponding to a doubled number of nodes along the 
Fermi surface. From the pairing amplitudes, we likewise 
observe that the pairing spreads out to the second nearest 
neighbor of the same sublattice. This is a consequence 
of the long-range Hubbard interactions: the Cooper pair 
wave function seeks to develop more nodes to minimize 
Coulomb repulsion, and is able to do so by longer range 
Cooper pairing. This, however, still does not tell us about 
the gap function of the d-wave instability. As the degen- 
eracy is protected by symmetry, the system could generi- 
cally form any linear combination dx2-y2 +e^^dxy of both 
d-wave solutions. For this purpose, we perform a mean 
field decoupling in the SC pairing channel and minimize 
the free energy as a function of the superposition param- 
eter. The necessary condition for such a minimum can 
be equivalently rephrased by satisfying the self-consistent 
gap equation [55] 
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The gap functions are displayed in Figs.jsj^aS) andjsjjbS). 
We always find d -I- id to be the energetically preferred 
combination. This is rather generic in a situation of de- 
generate nodal SC order parameters, since such a com- 
bination allows the system to avoid nodes in the gap 
function. The gaps we find are hence nodeless and only 
slightly change their anisotropy as the pairing function 
varies [Fig. [3^ and[3]3]. As graphene can be tuned rather 
accurately to the van Hove filling where we find the high- 
est critical scale, it may be a reasonably accesible exper- 
imental system to study such a SC phase. The expected 
experimental evidence for d + id would hence be a node- 
less gap detectable through transport measurements and 
singlet pairing due to a Knight shift drop below T^,. A mi- 
nor caveat is given by the role of impurities which may 
spoil the symmetry between the two d-wave solutions, 
which could give rise to a nodal gap beyond sufficient 
impurity concentration |26j . 

f-wave phase. It is likewise interesting to investigate 
the triplet /-wave instability |27j which dominates for 
longer ranged Coulomb interaction [Fig. |4| . It obeys the 
one dimensional Bi or B2 representation, depending on 
the range of the Coulomb interaction. For Ui/Uq = 0.5, 
the form factor and pairing amplitudes are plotted in 
as well as for Ui/Uq = 0.5, U2/U0 = 0.3 in 
We again find that the Cooper pair distance 
increases, which manifests itself as a change of the form 
factor /[/si] = sm{^/3ky) — 2 sin(-\/3fcj,/2) cos(3fc^/2) in 

(a) to /[/sa] — sin(3A:a;) — 2 sm{5kx/ 2) cos{3 V^ky/ 2) in 

(b) . The gap function follows the absolute value of the 
form factor, showing a nodal gap, where the points of 
the nodes change with increasing Coulomb range. In the 
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FIG. 4. (Color online). Pairing amplitudes, form factors, 
and gap profiles for the /-wave phase as defined in Fig. [3] 
representative for larger fillings than the VHS for (a) Ui/Uo = 
0.5 to (b) Ui/Uo = 0.5, U2/U= = 0.3. The gap profile is nodal, 
the nodal points shift from (a) to (b). 



case of /-wave, the position of the nodes would hence 
indicate the Cooper pairing distance associated with the 
long range properties of the Coulomb interaction, and 
suggest experimental evidence of a nodal gap from trans- 
port and an invariant Knight shift due to triplet pairing. 
For filling smaller than the VHS, the Fermi surface is 
disconnected and it can happen that the nodes do not 
coincide with the Fermi surfaces. While probably very 
low in Tc, depending on whether Bi or B2 is preferred, 
this /-wave regime could be nodeless. 

Summary and outlook. In summary, we have provided 
a detailed analysis of the competing many-body phases of 
graphene at and around van Hove filling. We find that for 
realistic band structure parameters and interactions, the 
exotic nodeless singlet d -|- zd-wave SC phase appears to 
be preferred over an extended phase space regime around 
the VHS. Variations of the kinetic parameters and effec- 
tive interaction scales can drive a transition to a spin 
density wave phase at the van Hove point. In a broader 
vicinity of the VHS, reduced Coulomb screening changes 
the form of the d + id Cooper pair wave function, and 
in certain limits might favor a nodal triplet /-wave SC 
phase. 

The possibility of the time-reversal symmetry breaking 
d + id phase in graphene is very intriguing: it has been 
noted early on in the context of the cuprates that such a 
phase would have various interesting properties such as 
quantized edge currents [551 HH] ■ Furthermore, provided 
Rashba spin-orbit interaction is present, d+id phase sup- 
ports Majorana modes in the vortex cores obeying non- 
Abelian statistics [3^ . The tunability of the Rashba in- 
teraction in graphene [3T] may enable realization of the 
Majorana modes; owing to the two-dimensional nature 
of graphene and its remarkable tunability, their obser- 
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vation and manipulation should be easier than in other 
materials. 
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